Mizer is a tool that can be used to simulate a dynamic size spectrum in an marine ecosystem, subject to changes through time, such as fishing pressure. Multiple interacting species and fishing gears can be defined allowing for a range of fisheries management strategy scenarios, and their ecosystem impacts, to be tested.

Below we explore some examples of what can be done using a previously calibrated Mizer model.

First, we load the necessary packages.

Species Collapse and Recovery in a Heavily Fished Multispecies System

In this example we will read in a calibrated model for the North Sea, which consists of 12 interacting species and a background resource. All species are predators and prey since they feed according to size-based rules and they encounter each other.

More information on how this model was calibrated with historical data can be found in our "mizerHowTo" package [tutorial link HTOM3] and in this [blogpost].

You can take a look at how the modelled species biomasses change through time by running the following:

# species' biomass through time
plotlyBiomass(sim,start_time = 1950,end_time = 2018)
# look at size spectra in final 5 years
plotSpectra(sim,time_range = c(2015:2018),power=2)
## Warning in validParams(params): You need to upgrade your MizerParams object with
## `upgradeParams()`.

Examining changes in the community relative to an unfished state

To be able to examine any changes in the community relative to an unfished state we can re-set effort = 0 and calculate the steady state.

## Simulation run did not converge after 99 years. Value returned by the distance function was: 0.172517288975348

We can compare the current size spectra (averaged over 2015-2019) to the unfished size spectra to assess whether there is any evidence of a size-structured trophic cascade due to fishing.

Here we can see the effect of the reduction in large sized individuals of heavily fished on the other sizes and species in the model, relative to the uifished steady state.

The abundance of some (but not all) of the smaller to medium sizes of prey are a lot higher when their larger predators are removed (note the logarithmic scale).

Sprat looks to be much lower.

We can do the same with Biomass to see when, if any, of the species collapse. For simplicity, we use < 0.1 of B/B_unfished as a proxy for a reference point for population collapse. We can add other reference points to this kind of plot, for example a simple rule of thumb for B_msy could be 0.5*B_unfished.

Projecting future recovery scenarios

We may now wish to explore the potential recovery the larger species and sizes in the system. To do this we set up another scenario, where the model starts with the last time step of the fished scenario.

First let's take a look at what the model has used in the historical period (from ICES stock assessments).

effort<-melt(sim@effort)  %>% filter(time >= 1980)

plot_effort<-ggplot(effort) + 
geom_line(aes(x = time, y = value, colour = gear))+ 
theme_pubr() +
labs(y="Effort",x="") 

plot_effort
## Warning: Removed 1 row(s) containing missing values (geom_path).

The species' maximum fishing mortality rates (which we call for convenience effort) for Cod and Saithe have declined but are not as low as what would be in line with single-species "Fmsy" (a reference level of fishing deemed sustainable from single-species stock assessements) and the fishing mortality rate for Sprat and Dab is very high.

We can the values directly like this:

sim@effort["2018",]
##      Sprat    Sandeel     N.pout        Dab    Herring    Gurnard       Sole 
## 1.24000000 0.08454123 0.31100000 1.27800000 0.19100000 0.37737415 0.33500000 
##    Whiting     Plaice    Haddock     Saithe        Cod 
## 0.19800000 0.19200000 0.24700000 0.41900000 0.32000136

Let's start a new simulation that begins with the effort from 2018 and projects forward for 50 years. We will apply a linear reduction in effort for a selected species to a target value (here assumed for simplicity to be 0.2 with effort expressed in terms of the species' fishing mortality rate for fully selected sizes).

To do this need to work with the effort array (time x gear) to enable changes in effort through time, for this scenario. Here we have a 'gear' for each species, since effort in this case were annual single-species fishing mortality estimates.

select_species="Sprat"
proj_effort<-sim@effort
target<-0.2

# reach target by 10 years
proj_effort[1:10,select_species]<-seq(from = proj_effort[1,select_species], to = target, length = 10)

# then hold at target
proj_effort[11:51,select_species]<-target

# check it
proj_effort[,select_species]
##         1847         1848         1849         1850         1851         1852 
## 2.930762e-27 2.222222e-02 4.444444e-02 6.666667e-02 8.888889e-02 1.111111e-01 
##         1853         1854         1855         1856         1857         1858 
## 1.333333e-01 1.555556e-01 1.777778e-01 2.000000e-01 2.000000e-01 2.000000e-01 
##         1859         1860         1861         1862         1863         1864 
## 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 
##         1865         1866         1867         1868         1869         1870 
## 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 
##         1871         1872         1873         1874         1875         1876 
## 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 
##         1877         1878         1879         1880         1881         1882 
## 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 
##         1883         1884         1885         1886         1887         1888 
## 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 
##         1889         1890         1891         1892         1893         1894 
## 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 2.000000e-01 
##         1895         1896         1897         1898         1899         1900 
## 2.000000e-01 2.000000e-01 2.000000e-01 3.479284e-16 5.736370e-16 9.457675e-16 
##         1901         1902         1903         1904         1905         1906 
## 1.559307e-15 2.570863e-15 4.238636e-15 6.988329e-15 1.152181e-14 1.899625e-14 
##         1907         1908         1909         1910         1911         1912 
## 3.131952e-14 5.163716e-14 8.513528e-14 1.403643e-13 2.314217e-13 3.815498e-13 
##         1913         1914         1915         1916         1917         1918 
## 6.290693e-13 1.037160e-12 1.709988e-12 2.819293e-12 4.648229e-12 7.663634e-12 
##         1919         1920         1921         1922         1923         1924 
## 1.263520e-11 2.083192e-11 3.434602e-11 5.662702e-11 9.336217e-11 1.539282e-10 
##         1925         1926         1927         1928         1929         1930 
## 2.537847e-10 4.184202e-10 6.898583e-10 1.137384e-09 1.875229e-09 3.091730e-09 
##         1931         1932         1933         1934         1935         1936 
## 5.097402e-09 8.404195e-09 1.385617e-08 2.284497e-08 3.766499e-08 6.209906e-08 
##         1937         1938         1939         1940         1941         1942 
## 1.023840e-07 1.688027e-07 2.783087e-07 4.588533e-07 7.565211e-07 1.247292e-06 
##         1943         1944         1945         1946         1947         1948 
## 2.056436e-06 3.390486e-06 5.589959e-06 9.216262e-06 1.519499e-05 2.505213e-05 
##         1949         1950         1951         1952         1953         1954 
## 4.130354e-05 6.809680e-05 1.122693e-04 1.850919e-04 3.051405e-04 5.030252e-04 
##         1955         1956         1957         1958         1959         1960 
## 8.291680e-04 1.366577e-03 2.251773e-03 3.708935e-03 6.105207e-03 1.003928e-02 
##         1961         1962         1963         1964         1965         1966 
## 1.648041e-02 2.697931e-02 4.396835e-02 7.113881e-02 1.137873e-01 1.788044e-01 
##         1967         1968         1969         1970         1971         1972 
## 2.736383e-01 4.034121e-01 5.663110e-01 7.500000e-01 9.336890e-01 1.096588e+00 
##         1973         1974         1975         1976         1977         1978 
## 1.226362e+00 1.321196e+00 1.605000e+00 1.647000e+00 1.439000e+00 9.510000e-01 
##         1979         1980         1981         1982         1983         1984 
## 6.230000e-01 2.152000e+00 1.041000e+00 9.540000e-01 1.602000e+00 9.170000e-01 
##         1985         1986         1987         1988         1989         1990 
## 1.263000e+00 1.040000e+00 3.500000e-01 1.241000e+00 3.020000e-01 1.468000e+00 
##         1991         1992         1993         1994         1995         1996 
## 6.850000e-01 8.000000e-01 1.557000e+00 6.780000e-01 1.439000e+00 1.330000e+00 
##         1997         1998         1999         2000         2001         2002 
## 9.760000e-01 1.645000e+00 8.620000e-01 1.422000e+00 1.494000e+00 1.577000e+00 
##         2003         2004         2005         2006         2007         2008 
## 1.207000e+00 1.845000e+00 1.272000e+00 1.618000e+00 1.559000e+00 1.311000e+00 
##         2009         2010         2011         2012         2013         2014 
## 8.460000e-01 9.760000e-01 8.840000e-01 1.347000e+00 1.190000e+00 5.650000e-01 
##         2015         2016         2017         2018         2019 
## 1.351000e+00 2.204000e+00 1.397000e+00 1.240000e+00 1.015000e+00
# run the simulation forward, using the 2018 abundances as initial values

sim_scen<-project(params,initial_n = sim@n["2018",,], initial_n_pp = sim@n_pp["2018",],effort=proj_effort)

# plot change in biomass relative to 2018 values

plotBiomass(sim_scen)

B_current<-getBiomass(sim_scen)[1,]

Brel_scen<-melt(sweep(getBiomass(sim_scen),2, B_current,"/"))

ggplot(Brel_scen)+ geom_line(aes(x=time,y=value,color=sp)) + 
geom_hline(yintercept = 1, linetype=1, colour="grey", size=0.75) +
theme_pubr()

We can see that when we reduce fishing on Sprat it increases the biomass of this species but also affects the biomass of other species in the community.

Are any species collapsed still?

Brel_ref<-melt(sweep(getBiomass(sim_scen),2, B0,"/"))

Brel_ref <- Brel_scen  %>% mutate(collapsed = value < 0.1)

ggplot(Brel_ref, aes(x = time, y = value, fill = collapsed)) +
  geom_col(position = "identity") +
  geom_hline(yintercept = 0.1, linetype=1, colour="red", size=0.75) +
   geom_hline(yintercept = 0.5, linetype=1, colour="dark grey", size=0.75) +
  scale_fill_manual(values = alpha(c("blue", "red"), .3)) +
  facet_wrap(vars(sp))+
  theme_void()

Exploration options - break into smaller groups:

  1. Cut and paste the above code chunks and edit it to explore your own fishing scenarios.

  2. Alter some of of the model parameters to see how this affects the scenario outcomes.

  3. Using the unfished sim object examine the effects of adding a species to the model

Further Exploration (Optional)

Example 2: Set up an explore a time-varying fishing scenario

We could examine how this fish community would have responded if there had been a completly different historical development of fishing fleets over time.

Here we will examine how some commonly used community indicators of the ecosystem effects of fishing change in response to this scenario.

We can alter the fishing parameters using a function called gear_params() and by changing the effort input.

Let's take a look at the fishing parameters. Note that the gears in the above model were already setup to be species-specific.

# look at the gear set up:
gear_params(params)
##                  species    gear       sel_func  l25  l50 catchability
## Sprat, Sprat       Sprat   Sprat sigmoid_length  7.6  8.1            1
## Sandeel, Sandeel Sandeel Sandeel sigmoid_length  9.8 11.8            1
## N.pout, N.pout    N.pout  N.pout sigmoid_length  8.7 12.2            1
## Dab, Dab             Dab     Dab sigmoid_length 10.1 20.8            1
## Herring, Herring Herring Herring sigmoid_length 11.5 17.0            1
## Gurnard, Gurnard Gurnard Gurnard sigmoid_length 19.8 29.0            1
## Sole, Sole          Sole    Sole sigmoid_length 16.4 25.8            1
## Whiting, Whiting Whiting Whiting sigmoid_length 19.8 29.0            1
## Plaice, Plaice    Plaice  Plaice sigmoid_length 11.5 17.0            1
## Haddock, Haddock Haddock Haddock sigmoid_length 19.1 24.3            1
## Saithe, Saithe    Saithe  Saithe sigmoid_length 13.2 22.9            1
## Cod, Cod             Cod     Cod sigmoid_length 35.3 43.6            1
##                  initial_effort
## Sprat, Sprat         1.29533333
## Sandeel, Sandeel     0.06510547
## N.pout, N.pout       0.31380000
## Dab, Dab             0.97800000
## Herring, Herring     0.18150000
## Gurnard, Gurnard     0.46250569
## Sole, Sole           0.37383333
## Whiting, Whiting     0.24266667
## Plaice, Plaice       0.18483333
## Haddock, Haddock     0.30150000
## Saithe, Saithe       0.39300000
## Cod, Cod             0.26666749

We can group these species together according to the gears they are caught by. Let's imagine a big super trawler.

# allocate species to gear types
gear_params(params)$gear <- c("super_trawler")
## The catchability has been commented and therefore will not be updated from the gear parameters.

Note that catchability is set to 1. This is because the fishing "effort" was here assumed to be the fishing mortality rate of fully selected sizes (see here link to setFishing on mizer website).

The previous effort won't work with these new gears, as it is gear x time. We only have a single gear now, so this is easier to set up.

Example 3: Adding a variable environment

Now we can look at probability of collapse.